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Abstract 

An analytical benchmark is proposed for graphene and carbon nanotubes, that may serve to test whatsoever 
molecular dynamics code implemented with REBO potentials. By exploiting the benchmark, we checked 
results produced by LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) when adopting 
the second generation Brenner potential, we made evident that the code in its current implementation 
produces results which are offset from those of the benchmark by a significant amount, and provide evidence 
of the reason. 
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1 Introduction 

Molecular dynamics (MD) simulations are nowadays more and more popular in scientific applications, es¬ 
pecially in those fields of material science involving nanotechnology and advanced material design. On one 
side, there are advantages in the speed and accuracy of the simulations, with the model of the potential for 
atomic interactions being optimized to reproduce either experimental values or quantities extimated by first 
principles calculations (considered, as a matter of facts, just like experimental results). On the other side, 
it is more and more frequent to use commercial or open-access codes implementing off-the-shelf potential 
models, and use them as a black box, without having a precise feeling with the code itself. One of the most 
used simulator is LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator), able to implement 
several interatomic potentials. By using an analytical discrete mechanical model, we present a benchmark 
for the equilibrium problem of graphene and carbon nanotubes, which can be applied to any kind of REBO 
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(reactive empirical bond-order) potential. The analytical condition proposed produces results in complete 
agreement with First Principles, Density Functional Theory and Monte Carlo simulations. With the aid 
of this benchmark, we show that LAMMPS code, when implemented with the second generation Brenner 
potential, produces results which are offset from those of the benchmark by a significant amount, and provide 
evidence of the reason. The purpose of this letter it is not to just to supply a test for the LAMMPS code, 
rather, it is to provide a general tool for testing any MD code. 


2 An analytical discrete model for equilibrium configurations of 
FGSs and CNTs 

The benchmark solution we propose has been developed within the context of carbon macromolecules, such 
as Flat Graphene Strips (FGSs) or Carbon Nanotubes (CNTs). When regarded from the point of view of 
MD, such aggregates are modelled as sets of mass points, whose configuration is described by the Cartesian 
coordinates of each point with respect to a chosen reference frame; each point is then interacting with the 
others - at least with the closest ones - and the interaction is captured by a suitable empirical potential, 
whose shape and parameters are htted with a set of selected experiments and ab initio calculations. The last 
generation potentials usually take into account multi-particle interactions, up to the third nearest neighbor, 
which is indispensable to capture the mechanics of complex systems, such as carbon macromolecules. 

In order to provide an easy-to-visualize mechanical picture, the perspective we here adopt is not the one 
of MD, we consider instead the approach of Favata et a/. [ 5 ], where a discrete mechanical model is detailed 
for 2 D carbon allotropes. In this view, the configuration of a molecular aggregate is not identified by the 
coordinates of the mass points, but rather by a suitable finite list of order parameters. In particular, the 
conditions of natural equilibrium of the aggregate can be determined and expressed in terms of such list and 
independently of the choice of the REBO potential. As we will see, the prediction of such equations are in 
total agreement with First Principles, Density Functional Theory and Monte Carlo simulations; moreover, 
given their generality, they can be exploited to establish benchmark solutions. 

In order to understand the physical meaning of the conditions we propose, we summarize some of the 
results of Favata et aL|^. Starting with the geometry, with reference to Fig. let the axes 1 and 2 be 
respectively aligned with the armchair and zigzag directions, and let ni, n2 be the number of hexagonal cells 
counted along these axes. Let us consider now the representative cell AiBiA2B^A^B2Ai. We note that 
the sides AiBi and A^B^ are aligned with the axis 1; the common length of corresponding bonds will be 
denoted by a, and we will call a-type the corresponding bonds. We see that the other four sides have equal 
length b {b-type bonds). We pass to introduce the bond angles and, since we intend to consider interactions 
up to the third neighbor, the dihedral angles. As to the bond angles, we notice that they can be of a-type 
and 13 -type (e.g., respectively, A3B2A1 and B2A1B1; see Fig. [^. As to the dihedral angles, there are hve 
types (01,..., ©5), which can be identified with the help of the colored bond chains in Fig. In conclusion, 
to determine the deformed configuration of a representative hexagonal cell, no matter if that cell belongs to 
a FGS or to an achiral GNT, we need to determine the 9 -entry order-parameter substring: 

isub'-= . ( 1 ) 

The complete order-parameter string for the whole molecular aggregate can be obtained by sequential jux¬ 
taposition of substrings. 


Due to the geometric compatibility conditions induced by the built-in symmetry (see Favata et al. for 
details), only three of the nine kinematic variables determine the natural configuration, which are chosen to 
be a, b, and a. In particular, by distinguishing the armchair (superscript A) from the zigzag (superscript Z) 
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a 



- ©1 - 6*3 

(or B2A1B1A2) 


02 ^ 3 - 02 ^ 1 ^] 


- 04 A^B^A^B'^ 

- 05 - 82 ^- 81^2 

(or B'^AiB^A^) 


Figure 1 : Order parameters in a graphene sheet. 


case, the order-parameter substring are given by, respectively: 

=ia,b,a,P {a,ip^), 

(a> 0^(a, 20 ^(q!, g>^), 0, 0); 

—'Z 

^sub (^7 ^5 ^5 {plT )5 

,^2 (a> 0> 2 of (a, 0). 


The explicit form of the functions /3 ’ , 0 (^, 0 ^’'^ is given in Favata et al.^ In (|^, = tt/tii is the angle 

between the plane of AiBiB^ and the plane of B1A2B2, when an armchair CNT is considered, and ip^ = ^ 
the angle between the planes of A1B1A2 and A2B3A3, when a zigzag CNT up is considered. In case of a 
FGS, we have = 0 , /3 ’ = tt — a/ 2 , and 0 ^ = 0 ^’’^ = 0 . 

The equilibrium equations turn out to be the following ones: 


0-0 = 0, o-fc = 0 , 

Ta + 213,a Tp + 01,a 7l -|- 202 ,q Tl + 03,Q Ts T 04,a Ta , 


where cfa,cfb,Ta,Tp, and %, are the so-called nanostresses, work-conjugate to changes of, respectively, bond 
lengths, bond angles, and dihedral angles of each type considered. The form of the third of @ depends on 
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which of the two achiral CNTs is dealt with; more precisely, we have that 

+ 20 ^,„ 7 ^^ + 0 ^,„ 7 ^-^ = 0 , 

rf + 2 / 3 ^,„ r| + 20f T/ + 0fT/ = 0 . 

Due to their generality, the conditions ([^ may serve as a benchmark for any REBO potential. To express 
the equilibrium equations in terms of the Lagrangian coordinates a, b, and a, it is necessary to introduce the 
constitutive equations for the stress, which result from the assignment of an intermolecular potential. In the 
next section, we detail the formulas in the Brenner 2 nd generation REBO potential [ 5 ] which are needed to 
solve (§ in terms of the order parameters. 


3 REBO potentials 

In the Brenner 2 nd generation REBO potential, the binding energy molecular aggregate is 

written as a sum over nearest neighbors: 

i J<I 

the interatomic potential Vjj is given by the construct 

Vij = Vji{rij) + bjjVAirij ), ( 4 ) 

where the individual effects of the repulsion and attraction functions Vji{rjj) and VA^rjj), which model 
pair-wise interactions of atoms 7 and J depending on their distance r/j, are modulated by the bond-order 
function bjj. The repulsion and attraction functions have the following form: 

VA{r) = f{r)'£Bne-^-\ 

n—1 

Vr{p) = Ae~'^'^ , 

where f^{r) is a cutoff function limiting the range of covalent interactions to nearest neighbors, and where 
Q, A, Bn, 7, and 6 n, are parameters to be chosen fit to a material-specific dataset. The remaining ingredient 
in @ is the bond-order function: 

= 2 > 

where superscripts a and n refer to two types of bonds: the strong covalent cr-bonds, between atoms in 
the same plane, and the 7r-bonds , which are perpendicular to the plane of cr-bonds. We now describe the 
functions b'ff^ and bjj. 

The role of function bff^ is to account for the local coordination of, and the bond angles relative to, 
atoms 7 and J, respectively; its form is: 

X! ffk{'riK)Gicos9ijK)e^‘-^‘^-\- 
\ Kyi,.j 

+ Pi.j{Nf,N^) 

Here, for each fixed pair of indices ( 7 , J), (a) the cutoff function fff^{r) limits the interactions of atom I to 
those with its nearest neighbors; (b) is a string of parameters designed to prevent attraction in some 

specific situations; (c) function Pjj depends on Nf and , the numbers of C and H atoms that are nearest 
neighbors of atom 7 , and is meant to adjust the bond-order function according to the specific environment 
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of the C atoms in one or another molecule; (d) according to Brenner et al.[2], for solid-state carbon, the 
values of both the string and the function Pjj are taken null. Finally, function G{cos9jjx) modulates 
the contribution of each nearest neighbor in terms of the cosine of the angle between the IJ and IK bonds; 
its analytic form is given by three six-order polynomial splines in cos 9, each of them defined in an interval 
of the bond angle 9. The corresponding coefficients are determined by fitting each polynomial spline to the 
values of G(cos 9) at certain values of 9. 

Function is given a split representation: 


b 


TT 

IJ 


= n 


RC 

IJ 


+ b 


DH. 
IJ I 


the first addendum depends on whether the bond between atoms / and J has a radical character and on 
whether it is part of a conjugated system, the second depends on dihedral angles. Function bfj^ is given by 

X [ X! X! J - ^IJKL)f?KjlK)fjL{TjL) 

\k(^I,J)L{^LJ) 

where function T[j is a tricubic spline depending on Nj = Nf + Np, Nj, and a function of local 

conjugation, while Qijkl is the diedral angle between the planes of I, J, K and /, J, L. 

When the point of view described in Sect. is assumed, the expressions of the potentials have to be 
specialized and written in terms of the order parameters in the substrings Q . On introducing the potentials 
Va and 14 for the a- and &-type bonds, we have, respectively: 

Va{a,l3, 0i) = VR{a) + ba{/3, 0i) VA(a), 

14(6, a,/3,02,03,04) = 14?(6) + 6f,(a,;3,02,03,04)14(6) 

(see Favata et a/.[B] for details). 

Once this has been done, the nanostresses entering the balance equations (§ can be expressed in terms 
of the order parameters by means of the following constitutive relations: 

= V^,(a) + 6a(40i)n(a), 

CTfe = 1/^(6) -|-6b(Q;,/3, 02,03,04)1 ^a(^) ! 

Ta = bb,a{a, /3, 02,03, 04) 14(6) , 

t/ 3 = ^ (4,/3(/3,0i) 14(a) + 26b,/3(a, /3,02, 03,04) 14(6)) , 

Ti = i6a,ei(/3,0i) 14(a), 

”4 = 2 bb,e2 (o^, J 02,03,04) 14(6), 

% = 66,03(0!, /3, 02,03, 04) 14(6) , 

Ti = 65,0^ (o, /3, 02,03, 04) 14(6) . 



4 Analytical vs LAMMPS results 


The most direct outcomes of our solution are natural geometry and energy, which can be used to check 
the correctness of whatever MD code. The results obtained by solving equations (|^ with Brenner 2nd 
generation potential are in good agreement with First Principles, Density Functional Theory (DFT) and 
Diffusion Monte Carlo (DMC) simulations, as Tables 1 and show. 

As an application of the possibility of exploiting (31 as a benchmark, we present in Table ^ the radii of 
a number of CNTs, showing that standard LAMMPS code underestimates them. In Table |4 the values of 
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the cohesive energy from our solution and those obtained with the use of LAMMPS code when adopting 
the 2nd-generation Brenner potential are presented; it can be seen that there is a remarkable difference: 
the cohesive energy is highly overestimated and our benchmark makes evident that the code in its current 
implementation definitely produces results which are offset from those of the benchmark. The origin of the 
discrepancies can be found only by a close inspection of LAMMPS source code. In fact, although in Brenner 
et aL[2] it is indicated that the values of the function Pjj should be taken null for solid-state carbon, the 
code assigns the value 0.027603. This latter value is actually dictated in Table VIII of Stuart et aZ.[Tn] for 
AIREBO potentials, due to the additional terms included in this potential. Whenever a LAMMPS user 
wants to adopt REBO potentials, he needs to change the hard-wired number for the variable PCCf[2] [0] in 
“pair_airebo.cpp”; unfortunately, the LAMMPS manual does not provide any information on this issue, and 
most studies based on LAMMPS REBO calculations are likely to have underestimation or overestimation 
of mechanical and geometrical properties presented in our Tables. An example of the use of LAMMPS 
with 2nd generation Brenner potential is Zhang et al. [U. When the value assigned in Brenner et aL[5] is 
implemented, the LAMMPS code produces the same results as the benchmark solution, letting alone a tiny 
difference due to numerical effects, as the third column of Tables [3| an d undeniably makes evident. 

Starting from the geometry and the energy gathered by means of (1^ it is possible to obtain secondary 
quantities. In Table and the Young moduli and the Poisson coefhcients are reported: the standard 
LAMMPS code overestimates the former and underestimates the latter. Our results are in very good agree¬ 
ment with the literature (see e.g. Agrawal et al. [I])- The differences between our benchmark and the 
LAMMPS code with modified parameters are ascribable to numerical effects, more accentuated because 
Young modulus and Poisson coefficients are quantities not directly evaluated, but rather derived, and an 
increment of numerical error is foreseeable. 
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Table 1: Cohesive eneregy (eV/atom) 

our El-Barbery Shin 

benchmark et al. [S] et al. 2014 [3] 

(First Principles) (DMC) 
-7.3951 n YH64 


Table 2: Radii (nm) of small CNTs, comparison with literature 


n, m) 

our 

benchmark 

Machon 
et al. 2002 [7] 
(DFT) 

Cabria 

et al. 2003 [4] 
(DFT) 

Popov 

2004 [g 

(TB) 

(3,3) 

0.211 

0.210 

0.212 

0.212 

(4,4) 

0.277 

- 

- 

- 

(5,0) 

0.208 

0.204 

0.206 

0.205 


Budyka 
et al. 2005 [3] 
(DFT) 


0.277 
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{n, m) 

Table 

Our 

benchmark 

(nm) 

3: Radii 
LAMMPS 
(standard) 
(nm) 

LAMMPS 

(modified) 

(nm) 

(3,3) 

0.2111 

0.2079 

0.2110 

(4,4) 

0.2767 

0.2723 

0.2766 

(5,5) 

0.3431 

0.3371 

0.3404 

(6,6) 

0.4101 

0.4035 

0.4100 

(7,7) 

0.4773 

0.4697 

0.4773 

(8,8) 

0.5447 

0.5361 

0.5447 

(10,10) 

0.6798 

0.6690 

0.6798 

(12,12) 

0.8151 

0.8022 

0.8151 

(18,18) 

1.2216 

1.2020 

1.2215 

(25,25) 

1.6961 

1.6689 

1.6960 

(5,0) 

0.2078 

0.2046 

0.2076 

(6,0) 

0.2447 

0.2409 

0.2446 

(7,0) 

0.2823 

0.2778 

0.2821 

(8,0) 

0.3202 

0.3151 

0.3201 

(9,0) 

0.3584 

0.3527 

0.3583 

(10,0) 

0.3969 

0.3905 

0.3967 

(12,0) 

0.4741 

0.4665 

0.4739 

(15,0) 

0.5905 

0.5810 

0.5904 

(20,0) 

0.7853 

0.7274 

0.7852 

(30,0) 

1.1760 

1.1572 

1.1759 



MD Benchmark: Testing LAMMPSs 


9 



Table 4: Cohesive energy 



Our 

LAMMPS 

LAMMPS 


benchmark 

(standard) 

(modified) 

(n, to) 

(eV / atom) 

(eV / atom) 

(eV/atom) 

(3,3) 

-7.0137 

-7.3838 

-7.0137 

(4,4) 

-7.1695 

-7.5569 

-7.1695 

(5,5) 

-7.2463 

-7.6422 

-7.2462 

(6,6) 

-7.2898 

-7.6905 

-7.2896 

(7,7) 

-7.3167 

-7.7204 

-7.3166 

(8,8) 

-7.3346 

-7.7403 

-7.3345 

(10,10) 

-7.3560 

-7.7640 

-7.3558 

(12,12) 

-7.3678 

-7.7771 

-7.3676 

(18,18) 

-7.3829 

-7.7038 

-7.3827 

(25,25) 

-7.3887 

-7.8003 

-7.3886 

(5,0) 

-6.9758 

-7.3417 

-6.9759 

(6,0) 

-7.0969 

-7.4763 

-7.0969 

(7,0) 

-7.1715 

-7.5593 

-7.1715 

(8,0) 

-7.2212 

-7.6144 

-7.2212 

(9,0) 

-7.2560 

-7.6531 

-7.2560 

(10,0) 

-7.2814 

-7.6812 

-7.2813 

(12,0) 

-7.3151 

-7.7186 

-7.3149 

(15,0) 

-7.3432 

-7.7499 

-7.3431 

(20,0) 

-7.3656 

-7.7747 

-7.3655 

(30,0) 

-7.3819 

-7.7927 

-7.3818 

graphene 

-7.3951 

-7.8074 

-7.3950 
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Table 5: Young modulus 



Our 

LAMMPS 

LAMMPS 


benchmark 

(standard) 

(modified) 

(n, to) 

(GPa) 

(GPa) 

(GPa) 

(3,3) 

893.9167 

987.0102 

885.0631 

(4,4) 

851.4536 

944.4810 

840.1351 

(5,5) 

804.6053 

901.869 

799.7630 

(6,6) 

800.9298 

891.6427 

789.8102 

(7,7) 

793.7379 

881.2895 

778.2888 

(8,8) 

784.8784 

872.9931 

767.9165 

(10,10) 

768.6379 

856.9461 

756.7625 

(12,12) 

756.3044 

846.2911 

746.5046 

(18,18) 

735.9332 

831.2163 

732.5252 

(25,25) 

726.2650 

823.3865 

726.4968 

(5,0) 

948.0854 

1046.3569 

943.9120 

(6,0) 

973.0048 

1075.4912 

968.5679 

(7,0) 

976.5265 

1082.0265 

971.5102 

(8,0) 

969.7954 

1076.6910 

965.8188 

(9,0) 

958.1824 

1066.6410 

954.9396 

(10,0) 

944.5151 

1053.1830 

941.3135 

(12,0) 

916.0510 

1025.8253 

915.1339 

(15,0) 

877.5820 

986.5745 

877.9641 

(20,0) 

830.0108 

940.5973 

835.7993 

(30,0) 

779.7321 

890.5147 

789.0504 
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Table 6: Poisson coefficient 



Our 

LAMMPS 

LAMMPS 

(n, to) 

benchmark 

(standard) 

(modified) 

(3,3) 

0.1450 

0.1237 

0.1563 

(4,4) 

0.2311 

0.2078 

0.2388 

(5,5) 

0.2924 

0.25782 

0.2963 

(6,6) 

0.3061 

0.27936 

0.3115 

(7,7) 

0.3181 

0.2937 

0.3318 

(8,8) 

0.3292 

0.3020 

0.3458 

(10,10) 

0.3466 

0.3194 

0.3526 

(12,12) 

0.3588 

0.3306 

0.3626 

(18,18) 

0.3779 

0.3426 

0.3740 

(25,25) 

0.3867 

0.3507 

0.3790 

(5,0) 

0.0655 

0.0362 

0.0661 

(6,0) 

0.0867 

0.0600 

0.0840 

(7,0) 

0.1100 

0.0800 

0.1105 

(8,0) 

0.1328 

0.1045 

0.1306 

(9,0) 

0.1544 

0.1234 

0.1511 

(10,0) 

0.1742 

0.1424 

0.1737 

(12,0) 

0.1923 

0.1725 

0.2032 

(15,0) 

0.2493 

0.2138 

0.2446 

(20,0) 

0.2950 

0.2564 

0.2835 

(30,0) 

0.3406 

0.2991 

0.3261 



